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Although genome-wide association studies (GWAS) have proven 
powerful for comprehending the genetic architecture of complex traits, 
they are challenged by a high dimension of single-nucleotide polymor¬ 
phisms (SNPs) as predictors, the presence of complex environmen¬ 
tal factors, and longitudinal or functional natures of many complex 
traits or diseases. To address these challenges, we propose a high¬ 
dimensional varying-coefhcient model for incorporating functional as¬ 
pects of phenotypic traits into GWAS to formulate a so-called func¬ 
tional GWAS or /GWAS. The Bayesian group lasso and the asso¬ 
ciated MCMG algorithms are developed to identify significant SNPs 
and estimate how they affect longitudinal traits through time-varying 
genetic actions. The model is generalized to analyze the genetic con¬ 
trol of complex traits using subject-specihc sparse longitudinal data. 

The statistical properties of the new model are investigated through 
simulation studies. We use the new model to analyze a real GWAS 
data set from the Framingham Heart Study, leading to the identifica¬ 
tion of several significant SNPs associated with age-specific changes 
of body mass index. The /GWAS model, equipped with the Bayesian 
group lasso, will provide a useful tool for genetic and developmental 
analysis of complex traits or diseases. 

1. Introduction. Phenotypic traits of paramount importance to agricul¬ 
ture and human health are quantitatively inherited, involving an unknown 
(usually very high) number of genes and undergoing a series of developmental 
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pathways and events [Lynch and Walsh (1998); Wu and Lin (2006)]. These 
complexities have made the genetic analysis of quantitative traits one of the 
most difficult tasks in biological sciences. Recently emerging genome-wide 
association studies (GWAS) have provided a great promise to systematically 
characterize the genetic control of complex traits and have been increasingly 
instrumental for the identification of significant genetic variants that con¬ 
trol phenotypic variation [Shuldiner et al. (2009); Takeuchi et al. (2009); 
Teichert et al. (2009); Yang et al. (2010)]. In human genetics, these results 
have started to gain a growing body of novel findings with potential clini¬ 
cal relevance [Daly (2010)]. In plant and animal genetics, GWAS, with the 
advent of a continuously falling genotyping cost, have been considered more 
seriously than any time before [Filiault and Maloof (2012)]. Despite their 
powerful impact on genetic studies, however, GWAS also encounter tremen¬ 
dous challenges from statistical analysis and interpretation. 

First, GWAS usually genotype hundreds of thousands of single-nucleotide 
polymorphisms (SNPs) on thousands of subjects, leading to a number of 
SNPs strikingly larger than the sample size used. Thus, to analyze these 
SNPs, simple univariate linear regression has to be used for individual tests. 
However, this method ignores the effects of other SNPs while assessing one 
particular SNP, and is subjected to a severe adjustment issue for multiple 
comparisons. Moreover, in biology and biomedicine, a phenotypic trait can 
always be better described by a dynamic trajectory because the trait un¬ 
dergoes a developmental process [Wu and Lin (2006)]. For example, human 
body height growth is a process from infancy to adulthood; the genetic study 
of adult height only, as conducted in many current GWAS [Lettre (2011)], 
provides limited information about the developmental genetics of height and 
its relationship with physical and mental characteristics at various stages of 
growth. In clinical trials, longitudinal measures are one of the most common 
data types, including HIV dynamics, cancer growth and drug response to 
varying doses [Wang et al. (2009)]. In this article, we address these issues by 
developing novel statistical models and algorithms that can analyze multiple 
SNPs simultaneously and integrate the developmental mechanisms of trait 
formation into a general GWAS framework through mathematical functions. 
The extension of the models to tackle genotype-environment interactions us¬ 
ing GWAS is straightforward. 

In a linear regression model for GWAS where SNPs are predictors, mul¬ 
tiple regression breaks down when the number of predictors far exceeds 
the number of subjects. Alternatively, variable selection approaches could 
identify important genetic factors and enhance the predictive power of the 
final model. For example, in analyzing case-control cohorts, lasso regression 
[Tibshirani (1996)] and elastic-net regression [Zou and Hastie (2005)] were 
studied by Wu et al. (2009) and Gho et al. (2009), respectively. Li et al. 
(2012) and He and Lin (2011) fnrther proposed two-stage variable selection 
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approaches to identify disease susceptibility genes. These methods, however, 
are restricted to models with a single phenotypic measurement from each 
subject. 

For genetic studies of dynamic traits that are measured repeatedly at mul¬ 
tiple time points, Wu and Lin (2006) proposed a conceptual model called 
functional mapping by incorporating longitudinal and functional data anal¬ 
ysis into a genetic design. Depending on the availability of explicit mathe¬ 
matical equations to describe a biological process, functional mapping uses 
parametric, nonparametric or semiparametric approaches for modeling non¬ 
linear effects of genetic variants over time and further revealing a dynamic 
landscape of interplay between genes and developmental pattern. Das et al. 
(2011) implemented functional mapping into a GWAS setting, leading to the 
birth of a so-called functional GWAS or /GWAS model. The basic principle 
of functional mapping and /GWAS is to model and predict the temporal 
pattern of genetic effects on a particular trait or disease in a quantitative 
manner. Time-varying change of gene expression has been found to be a 
ubiquitous phenomenon because different metabolic pathways, regulated by 
genes directly or indirectly, are required for an organism to best adapt to 
developmental alteration. In a genetic study of body mass index (BMI) by 
linkage mapping, Gorlova et al. (2003) identified different BMI susceptibil¬ 
ity genes as well as different modes of inheritance triggered by these genes 
in children and adults. A common variant in the obesity-associated FTO 
gene, identihed by a genome-wide search, was observed to be reproducibly 
associated with BMI and obesity from childhood into old age, but displayed 
varying magnitudes of genetic effects between child and adult stages [Fraying 
et al. (2007)]. 

To increase its applicability in clinical genomics, /GWAS could further 
accommodate irregular longitudinal data measured at subject-specihc time 
points. But both functional mapping and /GWAS analyze SNPs individu¬ 
ally or pairwise, and are incapable of depicting a comprehensive picture of 
the genetic architecture of dynamic traits. The motivation of this article is 
to develop a variable selection model for /GWAS, with a focus on nonpara¬ 
metric modeling of temporal genetic effects of SNPs. Variable selection in a 
nonparametrical setting is equivalent to selecting a subset of predictors with 
nonzero functional coefficients. Lin and Zhang (2006) developed GOSSO for 
model selection in a smoothing spline ANOVA model, with the penalty term 
being the sum of component norms. Zhang and Lin (2006) further extended 
it to nonparametric regression in an exponential family. Wang, Li and Huang 
(2008) estimated time-varying effects using basis expansion and selected sig- 
nihcant predictors by imposing SGAD penalty functions on the L 2 -norm of 
these basis expansions. 

We propose a Bayesian group lasso approach for variable selections in 
nonparametric varying-coefficient models. Group lasso was first proposed by 
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Yuan and Lin (2006). They considered the problem of selecting important 
groups of independent variables in linear regression models and generalized 
lasso by encouraging sparsity at the group level. However, since the Hessian 
is not defined at the optimal solution, they did not provide variance esti¬ 
mates for the regression coefficients. Here, we express time-varying effects 
as a linear combination of Legendre polynomials, and in such a case, the 
selection of important predictors corresponds to the selection of groups of 
polynomials. We develop a Bayesian hierarchical model for group variable 
selection and estimate all parameters by MCMC algorithms. Our method 
provides not only point estimates but also interval estimates of all param¬ 
eters, and the traditional Bayesian lasso [Park and Casella (2008)] is its 
special case in which the response variable is univariate. 

In Section 2, we introduce the /GWAS model that connects genotypes 
and irregular longitudinal phenotypical data. Section 3 shows the Bayesian 
hierarchical representation for this nonparametric varying-coefficient model, 
where group lasso penalties are applied to individual functional coefficients. 
The posterior computations as well as the interpretation of the results are 
described in Section 4. In Section 5, the statistical properties of the model 
are investigated through simulation studies. Section 6 provides the appli¬ 
cation to a real GWAS example from the Framingham Heart Study that 
analyzes age-specific changes of genetic effects on body mass index (BMI). 
BMI is a heuristic measure of body weight based on a person’s weight and 
height, providing the most widely used diagnostic tool to identify whether 
individuals are underweighted, overweighted or obese, and, further, to ex¬ 
amine their risk of developing obesity-related diseases, such as hyperten¬ 
sion, type 2 diabetes and cardiovascular diseases [Frayling (2007)]. We use 
a nonparametric approach based on orthogonal polynomials to approximate 
age-specific change in BMI. The discussion about the new model is given in 
Section 7. 

2. The fGWAS model. The model for functional genome-wide associ¬ 
ation studies (/GWAS) is the integration of functional data analysis and 
genome-wide association studies. The primary goal of the /GWAS is to 
study the dynamic pattern of genetic actions and interactions triggered by 
significant SNPs throughout the entire genome. Beyond traditional GWAS, 
/GWAS targets phenotypic traits that are measured longitudinally at re¬ 
peated time points. Suppose in a genome-wide association study involving 
n subjects, a continuous longitudinal trait of interest is measured at ir¬ 
regularly spaced time points, which are not common to all subjects. Let 
Yi = ... ,yi{tiTi))'^ be the Tj-dimensional vector of measurements on 

subject i where tj = {tn,... ,tiTi)'^ is the corresponding vector of measure- 
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ment time points after standardization, y* can be described as 

+ 


\yiitiTi), 


y{tiTi)/ \ai{tiTi) 


OiqitiTi)/ \Xiq/ 


( 2 . 1 ) 
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/ di{tii) ■ 
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/ ei{tii) \ 
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\di{tiTi) ■ 

dp{tiTi) / 
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+ 



We introduce matrix notation for a succinct presentation. Let OL{tig) = 
... ,a„{ti£))'^ be the g-dimensional vector of covariate effects, X* = 
{Xii,...,Xiqy be the observed covariate vector for subject i, aitu) = 
... ,Op(fi£))^ and d{tii) = (diitu ),.. ■,dp{tit)Y' be the p-dimensional 
vectors of the additive and dominant effects of SNPs, respectively. Further¬ 
more, let = (^ii) ■ ■ ■) iip)^ and = (^ji,..., Cip)^ be the indicator vectors 
of the additive and dominant effects of SNPs for subject i. Thus, at time 
point tit, 

yi{tii) = + cx{tiifXi + 

( 2 . 2 ) 


where is the overall mean and ei{tii) is the residual error assumed to 

follow a N(0, iT^(fi£)) distribution. The jth elements of and Ci are defined 
as 







if the genotype of SNP j is AA, 
if the genotype of SNP j is Aa, 
if the genotype of SNP j is aa, 
if the genotype of SNP j is Aa, 
if the genotype of SNP j is A A or aa. 


In other words, aj{tit) represents the average effect of substituting one allele 
for the other, and dj{tit) represents how the average genotypic value of the 
heterozygote deviates from the mean of the homozygotes. 

In the /GWAS model, the effects of covariates and SNPs are assumed to be 
functions of time. Many methods of estimating time-varying coefficients of a 
linear model in a longitudinal data setting have been proposed and studied, 
including basis expansion methods, local polynomial kernel methods and 
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smoothing spline methods. Among these techniques, Legendre polynomials 
have been widely used by quantitative geneticists for modeling the growth 
curves [Lin and Wu (2006)], the programmed cell death (PCD) process [Cui 
et al. (2008)] or the genetic effects responsible for other traits [e.g., Suchocki 
and Szyda (2011); Yang and Xu (2007); Das et al. (2011)]. By approximating 
time-varying effects using Legendre polynomials, the expansion coefficients 
can be solved through regression. Moreover, the biological evidence or the 
prior belief about the time-dependency of genetic control can be integrated 
by just truncating the series. Motivated by these studies, we approximate 
the effect of the kth covariate by a Legendre polynomial of order v — 1: 

(2.3) {ak{tii),...,ak{tiTi)) = Ljrfc, /c = l,...,g. 


where = (r^o, • • • are the Legendre polynomial coefficients, and 


(2.4) 



/I ta ^(3^-1) 

Vl ^(34.-1) 


are Legendre polynomial functions. Similarly, other time-varying effects can 
be represented as 


T 

(2.5) (aj(tii),...,aj(tirJ) =Uihj, j = l,...,p, 

T 

( 2 . 6 ) =UiCj, j = l,...,p, 

(2.7) = Uim, 

where hj = (fejO; • • • j &j(i;-i))^ are the Legendre polynomial coefficients for 
the additive effect of the jth SNP, Cj = (cjo, • • •, Cj(.u_i))'^ are the Legen¬ 
dre polynomial coefficients for the dominant effect of the jth SNP, and 
m = (mo,... ,m^_i)^ are the Legendre polynomial coefficients for the over¬ 
all mean function. 

After introducing Legendre polynomials to approximate time-varying ef¬ 
fects of covariates and SNPs, the full model of /GWAS becomes 


Viitii) = njiui + (u^ri,..., u^rg)X* 

(2.8) -L (u^bi,..., ufihp)$i + {ujici ,... ,u^Cp)Ci -L ei{ta), 




Last, since measurements within each subject are possibly correlated with 
one another, we assume that e* = (ej(tii),..., ei{tiT^))'^ follows a multivariate 
normal distribution with zero mean and covariance matrix Sj. Both para¬ 
metric and nonparametric methods have been developed to model the struc¬ 
ture of covariance between longitudinal measurements [Ma, Casella and Wu 
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(2002); Zhao et al. (2005); Yap, Fan and Wu (2009)]. In particular, we em¬ 
ploy the first-order autoregressive [AR(1)] model to approximate the residual 
covariance matrix. This covariance structure allows different measurement 
time points for different subjects, and assumes a constant variance over time 
and an exponentially decaying correlation, plU 2 -Li|^Q < ^ < between two 
measurements. Moreover, the matrix determinant in the likelihood function 
can be easily computed. In our real data example, the variance of repeated 
measurements is stable over time. In longitudinal data sets with variance 
heteroscedasticity, however, a Transform-Both-Sides (TBS) technique [Wu 
et al. (2004)1 can be employed to satisfy the variance stationarity assumption 
in the AR(1) model. 

3. Bayesian hierarchical representation for gronp Lasso penalties. In 

high-dimensional regression problems, such as GWAS, a regularized ap¬ 
proach is preferred to identify predictors with nonzero effects and to achieve 
better out-of-sample predictive performance. When parameters that we would 
like to penalize are finite-dimensional, we may apply different penalty func¬ 
tions to them to perform variable selection. But when these parameters are 
nonparametric smooth functions, a traditional regularization procedure can¬ 
not be directly applied. In this situation, regularized estimation for selecting 
important predictors is equivalent to selecting functional coefficients that are 
not identically zero. 

Let ||bjj| be the L 2 norm of the vector hj. The time-varying additive 
effect of the jth SNP is identically zero if and only if ||bj|| = 0. Therefore, if 
we estimate additive effects by a Legendre polynomial of order u, and would 
like to identify significant additive effects via penalized methods, we could 
partition all parameters of additive effects (b^, ..., bj) into p groups of size 
V according to p SNPs, and encourage sparse solution at the group level or 
select a subset of groups with nonzero L 2 norms. That is, the group lasso 
minimizes the following penalized least square: 

(3-1) + A^||m| + A*^||cj||, 

i=i i=i 

where = (yf ,..., y^), = Ey'^ = (/i^,..., /x^) and A and A* are two 

regularization parameters. A and A* control the amount of shrinkage toward 
zero; the larger their values, the greater the amount of shrinkage. They 
should be adaptively determined from the data to minimize an estimate of 
expected prediction error. 

From a Bayesian perspective, the group lasso estimates can be interpreted 
as posterior mode estimates when the regression parameters have multivari¬ 
ate independent and identical Laplace priors. Therefore, when group lasso 
penalties are imposed on the Legendre coefficients of additive and dominant 


LI, WANG, LI AND WU 


effects, the conditional prior for bj is a multivariate Laplace distribution 
with the scale parameter 

(3.2) 7r(bj|cj^) = (uA^/(7^)’'^^exp(—(uA^/(7^)~^'^^||bj||), 

and the conditional multivariate Laplace prior for dominant effect Cj is 

(3.3) 7r(cj|iT^) = (uA*^/(T^)^^^exp(—(uA*^/(T^)~^^^||cj||). 

To ensure the derived conditional distribution of bj has a standard form, 
we rewrite the multivariate Laplace prior distribution as a scale mixture of 
a multivariate Normal distribution with a Gamma distribution, that is, 

M-Laplace(bj|0, (uA^/cr^)”^^^) 

oc (uAVf^^)’'^^exp(-(uAVo-^)^''^||bj||) 

2 


OC 


J MVN(bj|0,diag((T^Tj ,... ,a'^Tj)) Gamma^Tj 


2 'vxy 


where (vX'^/a^) ^/^) is the scale parameter of the multivariate Laplace dis¬ 
tribution, a v-hy-v diagonal matrix diag(cr^rj ,..., cr^Tj ) is the covariance 
matrix of the multivariate normal distribution with mean zero, is the 
shape parameter of the Gamma distribution, and is the scale parameter 
of the Gamma distribution. After integrating out r?, the conditional prior 
on bj has the desired form (3.2). Then, in a Bayesian hierarchical model, 
we can rewrite the multivariate Laplace priors on bj as 

bjlr? cr^ ~ MVN(0,diag((T^T? ... 


Tj^|A ~ Gamma 


fv + 1 2 A 

' 2 ’uA2 j■ 




Likewise, the multivariate-Laplacian prior on Cj can be replaced by 


Cj|rj ,(T^ ~ MVN(0,diag((T"^T, 


,2^*2 
j ’■ 




T2|A~Gamma('^,^y 


Then, given A and A*, we have the following hierarchical representation of 
the penalized regression model: 


y|m,rfc,bj,Cj,p,cr^ oc (27r) EI*T)/2| jEj 


-1/2 l„-l/2Er(yi-/^d^S-l(yi-Md 


m~AA,(0,Smo), 
rfc ~ AA,(0, Sro), 


k = l,...,q, 
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f72,A,A* >0, 


where A and A* are regularization parameters or group lasso parameters that 
control the shrinkage intensities in estimating genetic effects. We assign a 
conjugate multivariate normal prior to m when estimating the overall mean 
function. We also assign conjugate multivariate normal priors to the Legen¬ 
dre coefficients of covariates r^, A: = 1,..., g, because covariates in GWAS are 
usually low dimensional and are not the target of variable selection. We as¬ 
sume a Uniform prior on [—1,1] for p, the autoregressive parameter in the as¬ 
sumed AR(1) covariance matrix. Finally, since the data are usually sufficient 
to estimate cr, we can use a noninformative prior such as vr(cr^) = 1/cr^ for cr^. 

Traditionally, two group lasso parameters A and A* can be prespecified 
by cross-validation or generalized cross-validation. However, in the Bayesian 
group lasso setting, A and A* can be estimated along with other parameters 
by assigning appropriate hyperpriors to them. This procedure determines 
the amount of regularization from the data and avoids refitting the model. 
In particular, the following conjugate gamma priors are considered. 




and 


where a, b, a* and b* are small values so that the priors are essentially 
noninformative. With this specification, group lasso parameters can simply 
join the other parameters in the Gibbs sampler. 

4. Posterior computation and interpretation. We estimate the unknown 
parameters and hyperparameters by sampling from their conditional poste¬ 
rior distributions through MGMC algorithms. Given the data likelihood and 
prior distributions, the posterior distributions of all unknowns can be ob¬ 
tained by Bayes’ theorem. For most of the parameters, the conditional poste¬ 
rior distributions have closed forms by conjugacy, which facilitates drawing 
posterior samples. 
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Assuming that priors for different predictors are independent, we can 
express the joint posterior distribution of all parameters as 

7r(m, Ffc, bj, , A, cj , rf ,X*,a‘^,p\y) 

q 

oc 7r(y|-)7r(m)7r(cr^)7r(p) n 

k=l 

X YlTr{hj\Tf)7r{T^\X)7r{X)7r{cj\Tf)Tr{Tf\X*)Tr{X*). 

i=i 

Conditional on the parameters (r^,bj,r?,A,Cj,A*,cr^,p), we derive 
the conditional posterior distribution of m as 


7r(m|y, r^, bj, r|, A, A*, c^^ p) 

0C7r(m)7r(y|-) 

( ^ T^p-l 

ocexpl --m 


oc 


oc 


“ ^ - UiinfT.^ ^(yi - - Uim) j 

i=i J 

exp I 

V i=i 

- 2^(yi - J 

i=l ) 

/ / n \ n ^ 

exp ( m'^ I + '^Uf j m - 2 ^(y* - 




2=1 


2=1 


Hence, the conditional posterior distribution of m is MVN^(/x^, Sm), where 


2=1 


V. 2=1 


and 




2 = 1 


-1 


BAYESIAN GROUP LASSO AND FUNCTIONAL GWAS 


11 


Similarly, since r^, hj and Cj have conjugate multivariate normal priors, 
the posterior distribution for is MVN„(/^^^,with 


2 = 1 




^ 2=1 


and 




the posterior distribution for bj is MVN„(/i(,^, ), with 


-1 




V. 2=1 


and 


-1 




2=1 


and the posterior distribution for Cj is MVN„(//^^., Sc^.), with 


-1 


Me, = 


X fE(yi-M*(-c,)fsri(C,,C/0^ , 


^ 2 = 1 


and 


-1 


Sc, = (uVf)-' +5;(Cut^*)^Sri(C,,[/0 


2=1 


Now, we derive the conditional posterior distribution for tJ and from 
the joint posterior distribution. Since 


7r{T^\y, m, , b,-, A, c,-, rf ,X*,a‘^,p) 
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and 


oc Tr{Tj\X)Tr{hj\Tj ,a^) 




expl 


xexp( -Xbj{a^diag{Tj,...,Tf)) ^bj 


cx exp ( —Tj 


;A2 


2aVf' ^ 




Tr(X |y,m,rfc,bj,T-,Cj,r* ,A*,f7 ,p) 


oc TT 


i=i 


^ /7A2\ (’^+l)/2 

oc (A^)“ ^exp(—6A^)^ 
j=i 


2 7 


z;A^ 2 
exp I -—r, 


the posterior distribution for is inverse-Gaussian(uA^, and the 


,-2 

posterior distribution for A^ is Gamma(a + 21^, 6-|- 

Similarly, the posterior distribution for is inverse-Gaussian(uA*^, 

b 

and the posterior distribution for A*^ is Gamma(a* + , b* + 

vJ2^- 

— ^ 2 ^ ^ )■ From these posteriors, we can see that the hierarchical expan¬ 
sion of the Multivariate Laplace prior indeed gives closed forms of posterior 
distributions for efficient Gibbs sampling. 

Last, if we assume a stationary AR(1) covariance structure, that is, 

/ 1 p\Ul-ti2\ ... p\tiTi-Ul\\ 

■)\ti2-Ul\ I ... p\UTi-U2\ 


Si = cr^Fi = cr^ 


YpiLi-L tJ ^|L2-LtJ 


1 


/ 


the posterior distribution for is an inverse chi-square distribution, or 


7r(cj2|-) ~Inv-xM 


^ i=l 


Er=i(y»-Ri)^r. \yi-tii 


where the first parameter is the degree of the freedom parameter and the 
second one is the scale parameter, and 

7 r(p|-) 0 C 7 r(y|-) 7 r(p) 
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2=1 


n 



Based on this expression, the corresponding Metropolis-Hastings algorithm 
can be developed to update p. 

We use MCMC algorithms to estimate the posterior distribution of each 
parameter by drawing posterior samples from the corresponding conditional 
posterior distribution, given the current values of all other parameters and 
the observed data. We use the potential scale reduction factor [PSRF; Gel- 
man and Rubin (1992); Gelman et al. (2004)] to access the convergence. 
Squared PSRF is defined as the ratio of the marginal posterior variance to 
the within-chain variance, and a PSRF less than 1.1 indicates good conver¬ 
gence. We run 4000 additional iterations after all chains converge. 

5. Computer simulation. We hrst investigate the new Bayesian group 
lasso approach for selecting important time-varying effects through simula¬ 
tion studies. We generate data in the /GWAS setting according to the model 
(2.8) with the number of covariates q = 1, the number of SNPs p = 3000, 
and the number of individuals n = 600 or 800. Following the simulation 
techniques in the literature, genotypical data ^ij is derived from Uij for 
i = 1,..., n and j = 1, ... ,p, where each Uij has a standard normal distribu¬ 
tion marginally, and cov{uij,Uik) = pc = 0.1 or 0.5, representing two levels 
of linkage disequilibrium. We set 


1, Uij ^ C, 

0, —c < Uij < c. 

1, Uij C, 



'ij ^ 


where c is used to determine the minor allele frequencies. Then, we derive 
the indicator matrix Qj of dominant effects from ^ij. 

We assume that the dynamic pattern of the trait is controlled by 5 SNPs 
and 1 covariate. In particular, we set hj = 0 for / = 4,... ,p, and Cj = 0 for 
/ = 1,2,6,... ,p. Sex is included as a covariate and is generated by randomly 
assigning a sex to each subject. The time-varying effects of overall mean, 
covariate and causal SNPs are generated by Legendre polynomials, with 
Legendre coefficients listed in Table 1. The true polynomial degrees for these 
causal SNPs could be 0, 1, 2 or 3, allowing constant genetic effects, linear 
genetic effects or more complicated patterns of genetic control. 

To simulate irregular longitudinal phenotypical data, we assume that the 
number of measurements for each subject is between 5 and 12, and all sub¬ 
jects are in the age range of 30 to 80 years. For each subject with a specihc 
number of measurements, traits of interest are observed at ages randomly 
drawn from 30 to 80. The residual covariance matrix among different time 
points was assumed to be AR(1) with p = 0.4 and cr^ = 4,9 or 16. The pheno- 
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Table 1 

Parameters used in the simulated example 


Time-varying effect 

Parameter 

0 

Legendre coefficients 

1 2 

3 

Mean effect 

m 

13.40 

-3.08 

1.88 

-3.20 

Covariate effect 

ri 

3.00 

0.15 

-2.67 

3.25 

Additive effect 

bi 

1.04 

0.88 

-2.05 

0.00 


b2 

1.17 

-0.22 

0.74 

-4.72 


ba 

1.40 

0.00 

0.00 

0.00 

Dominant effect 

C 3 

1.49 

-2.13 

4.82 

1.42 


C4 

1.00 

1.32 

1.90 

1.50 


C5 

1.26 

-1.22 

0.00 

0.00 


types observed at subject-specific time points and genotypes of all subjects 
are collected for Bayesian analysis. 

For each simulated data set, we implement MCMC algorithms as de¬ 
scribed in Section 4. In practice, the degree of Legendre polynomials should 
be determined a priori. We recommend a procedure that analyzes all SNPs 
with different polynomial degrees, where group lasso penalties are used to 
regularize the estimation. When the polynomial degree is 0 (constant effect), 
the group lasso penalty reduces to a lasso penalty. Then the polynomial de¬ 
gree V that gives the lowest Bayesian information criterion (BIC) of the final 
model is chosen. In simulations, however, this is computationally expensive. 
Therefore, the polynomial degree is fixed at h = 3 in simulation studies. Sim¬ 
ulation results (see Table 2) suggest that, as long as the specified polynomial 
degree is greater than or equal to the largest degree of all nonzero effects, 
the proposed framework works well in selecting casual SNPs and estimating 
their time-varying effects. 

Once all posterior samples are collected from MCMC algorithms, SNPs 
are selected in the following way: a time-varying additive effect aj{t) or 
dominant effect dj[t) is included in the final model if at least one of its four 
Legendre coefficients has a two-sided 95% interval estimate that does not 
cover zero. In the supplemental article [Li et al. (2015)], we plot the potential 
scale reduction factor against iterations for each parameter in bi, b 2 , bs, 
C3, C4 and C5. This is a simulation randomly drawn from the specification 
n = 600 and = 16. All chains converge very quickly and stay below the 
threshold of 1.05 (the red line). 

To evaluate the variable selection performance of the proposed procedure, 
we calculate several measures of model sparsity for the final model, which 
are summarized in Table 2. Column “C” shows the average number of SNPs 
with nonzero varying-coefficients correctly included in the final model, and 
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Table 2 

Variable selection performance in the simulated example 




No. of nonzeros 


Proportion of 



n 


C 

IC 

Under-fit 

Correct-fit 

Over-fit 

Time (h) 

pa = 
600 

0.1 

16 

3.77 

0.00 

0.86 

0.14 

0.00 

17.99 

600 

9 

4.93 

0.00 

0.07 

0.93 

0.00 

17.67 

600 

4 

5.00 

0.00 

0.00 

1.00 

0.00 

17.40 

800 

16 

4.99 

0.00 

0.01 

0.99 

0.00 

23.69 

800 

9 

5.00 

0.00 

0.00 

1.00 

0.00 

24.78 

800 

4 

5.00 

0.00 

0.00 

1.00 

0.00 

24.35 

pa = 
600 

0.5 

16 

4.61 

0.00 

0.35 

0.65 

0.00 

17.90 

600 

9 

5.00 

0.00 

0.00 

1.00 

0.00 

17.29 

600 

4 

5.00 

0.00 

0.00 

1.00 

0.00 

17.63 

800 

16 

5.00 

0.00 

0.00 

1.00 

0.00 

23.97 

800 

9 

5.00 

0.00 

0.00 

1.00 

0.00 

23.49 

800 

4 

5.00 

0.00 

0.00 

1.00 

0.00 

24.38 


column “IC” is the average number of SNPs with no genetic effect incorrectly 
included in the final model. Column “Under-fit” represents the proportion of 
excluding any relevant SNP in the final model. Similarly, column “Correct- 
fit” represents the proportion that the extract true model was selected and 
column “Over-fit” gives the proportion of including all relevant SNPs as 
well as one or more irrelevant SNPs. Clearly, both sample size and the noise 
level play important roles in how well the Bayesian group lasso could se¬ 
lect the exactly correct model. However, as sample size decreases and noise 
increases, our procedure tends to select fewer important SNPs rather than 
produce more false positives. Moreover, the impact of linkage disequilibrium 
is limited, and our method works slightly better in the presence of high 
linkage disequilibrium. 

Other than the performance of selecting truly important SNPs, we fur¬ 
ther investigate how well the procedure estimates the time-varying effects 
of selected SNPs. To ameliorate the bias of the parameter estimates intro¬ 
duced by group lasso penalties, we always refit the /GWAS model after 
variable selections, where only selected SNPs are included in the final model 
and all regularization parameters are set to zero. For each time-varying ge¬ 
netic effect of important SNPs, Tables 1 and 2 in the supplemental article 
[Li et al. (2015)] summarize the average estimates, standard errors and the 
mean squared errors (MSEs) of Legendre coefficients over replications where 
the effect is selected for pc = 0.1. As can be seen from these tables, both 
bias and standard error decrease as noise level decreases. MSEs are slightly 
lower for additive effects and lower order Legendre coefficients. 
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To compare the parameter estimates with those produced by another 
strategy aimed at the same genetic model, we implement the univariate 
/GWAS approach by Das et ah (2011) using the same data set. Specifically, 
this single-SNP analysis extends the traditional GWAS analysis framework 
by allowing the phenotype to be collected repeatedly over time and ap¬ 
proximating the time-varying genetic effects by Legendre polynomials. A 
Benjamini-Hochberg false discovery rate (FDR) controlling procedure is used 
to adjust for multiple comparisons in selecting signihcant SNPs. Table 3 in 
the supplemental article [Li et al. (2015)] shows that this single-SNP analysis 
produces biased estimates for all parameters.^ 

Finally, we compare the variable selection performance of four approaches: 
(1) a Bayesian group lasso; (2) a univariate /GWAS approach by Das et al. 
(2011); (3) afunctional principal component analysis (fPGA) approach [Ram¬ 
say and Silverman (2005)] that analyzes the fPCA of the longitudinal phe¬ 
notype; and (4) a slope model that simplifies the longitudinal phenotype to 
its slope. ^ In the third and the fourth model, the leading three fPCA scores 
and the slope calculated from each growth curve are tested against genetic 
predictors, respectively, where group lasso or lasso regressions with 5-fold 
cross-validation are used to select relevant SNPs. 

For fairness of comparison, longitudinal phenotype data are not gen¬ 
erated from our nonparametric genetic model (2.8). Instead, we use the 
same genotype data with pc = 0.1 but assume the following time-varying 
genetic effects: ai{t) = 0.5 -|- sin(0.2t), a 2 {t) = 1/(0.5 -|- exp(—O.OGt)) — 0.5, 
a 3 (t) = log(0.05t), dsit) = —1.5, d 4 ^{t) = 60/t, and d 5 {t) = 0.2 — 0.035t for 
the first five SNPs. These functional forms are unknown to researchers. Ta¬ 
ble 3 presents variable selection results, where all measures strongly prefer 
the Bayesian group lasso. Among the alternative approaches, the univariate 
/GWAS approach has the best variable selection performance. For the fPCA 
approach and the slope approach, the probability of selecting casual SNPs 
increases with the signal-to-noise ratio (column “C”), but the proportion 
of under-fit is always substantial. Interestingly, as signal-to-noise ratio in¬ 
creases, the probability of identifying false positives also increases (columns 
“IC” and “Over-fit”), especially when decreases from 16 to 9. The incon¬ 
sistency of these procedures suggests the risk of inflated false positive rates 
when only the major movements of growth curves are captured and tested 
in association studies. 

In the above simulation studies, the minor allele frequency is set to 0.3. 
Unreported simulations also demonstrate that as the minor allele frequency 


^ Since this approach cannot identify if the significance is due to the additive effect or 
the dominant effect, both effects are reported for five important SNPs. 

®We thank the Associate Editor and an anonymous referee for pointing out the fPCA 
method and the slope method, respectively. 
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Table 3 

Variable selection performance of alternative methods in the simulated example 




Nonzeros 

Proportion of 

Nonzeros 

Proportion of 

n 

rr" 

C 

IC 

U.-fit 

C.-fit 

O.-fit 

C 

IC 

U.-fit 

C.-fit 

O.-fit 

Bayesi 

600 

ian 

16 

group lasso 

3.93 0.00 

0.83 

0.17 

0.00 

Das et al. (2011) 

4.94 0.51 0.06 

0.55 

0.39 

600 

9 

4.80 

0.00 

0.20 

0.80 

0.00 

5.00 

0.20 

0.00 

0.83 

0.17 

600 

4 

5.00 

0.00 

0.00 

1.00 

0.00 

4.98 

0.02 

0.02 

0.96 

0.02 

800 

16 

4.86 

0.00 

0.14 

0.86 

0.00 

4.99 

0.40 

0.01 

0.70 

0.29 

800 

9 

5.00 

0.00 

0.00 

1.00 

0.00 

5.00 

0.23 

0.00 

0.81 

0.19 

800 

4 

5.00 

0.00 

0.00 

1.00 

0.00 

5.00 

0.01 

0.00 

0.99 

0.01 

Functional PCA 
600 16 0.51 

4.46 

1.00 

0.00 

0.00 

Slope model 
1.06 10.28 

1.00 

0.00 

0.00 

600 

9 

1.19 

7.01 

0.98 

0.00 

0.02 

2.53 

16.42 

0.99 

0.00 

0.01 

600 

4 

2.75 

13.53 

0.78 

0.00 

0.22 

3.65 

21.22 

0.97 

0.00 

0.03 

800 

16 

0.77 

5.02 

1.00 

0.00 

0.00 

1.82 

10.91 

1.00 

0.00 

0.00 

800 

9 

1.79 

11.14 

0.88 

0.00 

0.12 

3.06 

14.90 

0.98 

0.00 

0.02 

800 

4 

2.90 

17.16 

0.61 

0.00 

0.39 

3.82 

19.00 

0.99 

0.00 

0.01 


decreases, both statistical powers and false positive rates decrease. But 
our method is still much better than the alternative approaches. Despite 
the Bayesian framework’s theoretical advantages in handling parameter un¬ 
certainty, practically it could be slower than frequentist methods. When 
n = 600, = 9, pg = 0.1 and the number of SNPs p = 1000, the Gibbs 

sampler’s computational time is about 5.70 hours. Experiments show that a 
linear regression line® can describe almost perfectly the relationship between 
the computational time in hours and p\ \ogiQ{time) = 0.754-|-logio(p/1000). 

6. Worked example. We use the newly developed model to analyze a real 
GWAS data set from the Framingham Heart Study (FHS), a cardiovascu¬ 
lar study based in Framingham, Massachusetts, supported by the National 
Heart, Lung, and Blood Institute, in collaboration with Boston University 
[Dawber, Meadors and Moore (1951)]. Recently, 550,000 SNPs have been 
genotyped for the entire Framingham cohort [Jaquish (2007)], from which 
493 males and 372 females were randomly chosen for our data analysis. These 
subjects were measured for body mass index (BMI) at multiple time points 
from age 29 to age 61. The number of measurements for a subject ranges 
from 2 to 18, and the intervals of measurement are also highly variable 
among subjects. As is standard practice, SNPs with rare allele frequency 
<10% were excluded from data analysis. The numbers and percentages of 


®We thank the Editor for sharing the idea of using this regression. 
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nonrare allele SNPs vary among different chromosomes and range from 4417 
to 28,771 and from 0.64 to 0.72, respectively. 

A single-SNP analysis was used to analyze the phenotypic data of BMI 
for males and females separately. Figure 1 gives — log^Q p-values for each 
SNP in the two sexes, from which 33,239 SNPs with — log;^o p-values greater 
than 2.0 in at least one sex were selected. Before applying Bayesian group 
lasso analysis to this irregular longitudinal data set, we imputed missing 
genotypes for a small proportion of SNPs according to the distribution of 
genotypes in the population. Then, by treating the sex as a covariate, we 
imposed group lasso penalties on both additive effects and dominant effects 
in hopes of identifying SNPs with notable effects on BMI, where all effects 
are possibly functions of time. According to our discussions in Section 5, 


Mate 



1 2 3 4 5 6 7 e 9 10 11 12 13 14 16 IS 17 1 8 20 22 


Female 



1 2 3 4 6 6 7 0 9 10 11 12 13 14 15 10 17 1 0 20 22 

Chfomoscmi* 


Fig. 1. Manhattan plot of p-values for association by genomic position for male and 
female, where different colors across the x-axis represent different chromosomes, and the 
horizontal line indicates the significance level obtained by the Benjamini-Hochberg FDR 
adjustment at a = 5%. 
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the whole procedure was repeated with polynomial degrees: 0, 1, 2, 3 and 4, 
and the corresponding BICs of the final model are as follows: 27,470, 27,444, 
27,416, 27,408 and 27,426. Therefore, a polynomial degree of 3 is appropriate 
in this real data example. 

The Bayesian group lasso selected 24 significant SNPs, located on chro¬ 
mosomes 1, 2, 3, 4, 6, 7, 12, 14, 16 and 23. Table 4 tabulates the names, 
positions, alleles and estimated Legendre coefficients of these SNPs. The 
first allele in the column “Alleles” represents the minor allele. Using the 
Legendre coefficient estimates, we plot their time-varying additive effects 
and dominant effects in Figures 2 and 3, respectively, where the associated 
interval estimates^ are also provided. Some of these detected SNPs are lo¬ 
cated in a similar region of candidate genes for obesity. For example, the 
detected SNPs on chromosomes 4, 6 and 12 are close to candidate genes for 
BMI-related type 2 diabetes [Frayling (2007)]. 

Figures 2 and 3 show that the time courses of the genetic effects of some 
SNPs are relatively constant (magenta), monotonically increasing (black) or 
decreasing (blue). That is, given a population carrying one of these SNPs 
in the same environment, the expected BMI is different at different ages. 
Individuals carrying certain SNPs may have lower BMI in mid-life but tend 
to have higher BMI when they are younger or older (red). Conversely, indi¬ 
viduals carrying certain SNPs tend to have higher BMI in mid-life (green), 
which may increase the risk for stroke later in life, according to a prospective 
study [Jood et al. (2004)]. 

7 . Discussion. When the number of predictors p is much larger than the 
number of observations n, highly regularized approaches, such as penalized 
regression models, are favorable to identify nonzero coefficients, to enhance 
model predictability and to avoid overfitting [Hastie, Tibshirani and Fried¬ 
man (2009)]. In this article, we proposed a Bayesian regularized estimation 
procedure for nonparametric varying-coefficient models that could simulta¬ 
neously estimate time-varying effects and implement variable selection. The 
procedure extends the standard Bayesian lasso [Park and Casella (2008)] 
and standard group lasso [Yuan and Lin (2006)] to a nonparametric setting, 
and is applicable to irregular longitudinal data. 

We approximated time-varying effects by Legendre polynomials and pre¬ 
sented a Bayesian hierarchical model with group lasso penalties that en¬ 
courages sparse solutions at the group level. The group lasso penalties are 


^Suppose for one varying-coefficient, the interval estimate of the gth Legendre coef¬ 
ficient is (bq,( 7 , </= 1, ■ ■ ■,4, and the Legendre polynomials are (mo,mi,M 2 ,M 3 )^ = 

(l,t, — 1), |(5U — 3t))^ for each standardized time point t € [—1,1]. Then the in¬ 
terval estimate of the varying-coefficient at time t is 69 ,( 71 ^ 9 , Y]q=i ^ 9 ,where 

bq,u = 69,(7 if Uq is positive and bq^L otherwise, and bq^L ~ bq^L if Uq is positive and bq^u 
otherwise. 



to 

o 


Table 4 

Information about selected SNPs in the real data example 

Estimated legendre coefficients 


Chr 

Name 

Position 

Alleles 


Additive effect 



Dominant effect 


1 

SS66334458 

79,393,823 

C/T 

-1.504 

-1.656 

-1.334 

-0.143 

1.136 

2.645 

2.480 

0.780 

1 

SS66050888 

93,240,623 

A/G 

-0.431 

-0.532 

-0.344 

-0.111 

0.737 

0.745 

-0.185 

0.190 

1 

SS66275851 

93,245,738 

C/T 

-0.128 

-0.949 

-1.096 

-0.448 

0.079 

-0.254 

0.278 

-0.004 

1 

SS66048018 

115,427,398 

A/G 

0.396 

0.217 

0.028 

0.418 

0.213 

-0.007 

0.571 

0.294 

1 

SS66287256 

221,051,934 

G/A 

0.497 

0.788 

0.934 

-0.047 

0.386 

-1.065 

-0.672 

0.221 

1 

SS66104828 

234,701,498 

A/C 

0.111 

-0.620 

-0.951 

-0.461 

1.445 

1.833 

1.307 

0.205 

2 

SS66484730 

103,489,666 

G/A 

-0.341 

0.254 

0.335 

0.098 

0.057 

0.612 

0.552 

-0.565 

2 

SS66232775 

103,493,541 

T/C 

0.476 

-1.098 

-0.806 

-0.220 

0.043 

-0.816 

-1.011 

0.810 

2 

SS66185516 

239,065,169 

G/T 

0.397 

0.074 

-0.053 

0.228 

1.039 

0.852 

0.687 

0.269 

3 

SS66397464 

73,251,862 

C/T 

0.415 

0.183 

-0.198 

-0.192 

0.677 

0.895 

0.437 

-0.212 

4 

SS66402098 

186,281,132 

T/C 

-0.225 

0.630 

0.565 

-0.009 

0.418 

0.651 

0.744 

0.244 

6 

SS66218814 

3,311,818 

C/T 

-0.724 

0.043 

0.159 

0.182 

0.237 

-0.795 

-1.056 

-0.336 

7 

SS66083459 

89,430,534 

T/G 

-0.744 

0.518 

0.336 

0.141 

0.377 

-1.070 

-1.555 

-1.214 

12 

SS66288005 

29,860,263 

A/G 

-0.342 

0.724 

0.541 

0.322 

0.096 

0.563 

0.875 

0.806 

14 

SS66282595 

24,339,998 

G/A 

1.461 

1.471 

1.217 

-0.246 

-0.588 

-1.194 

-0.973 

0.022 

14 

SS66411959 

24,340,175 

G/A 

-0.782 

-1.357 

-1.311 

-0.080 

0.307 

0.323 

0.068 

-0.276 

14 

SS66416767 

24,348,496 

G/T 

-0.232 

-0.589 

-0.117 

-0.033 

0.507 

0.610 

-0.098 

-0.559 

14 

SS66281419 

77,702,561 

G/A 

-0.802 

0.151 

0.488 

0.252 

0.438 

-0.109 

0.254 

-0.189 

16 

SS66091573 

57,829,089 

C/T 

1.402 

2.674 

1.450 

0.400 

0.999 

2.747 

1.859 

0.426 

16 

SS66242525 

57,935,351 

C/T 

-0.548 

-0.516 

0.465 

0.579 

-0.537 

-0.479 

0.093 

1.234 

16 

SS66489647 

57,938,934 

A/G 

-0.217 

-2.058 

-1.834 

-1.059 

0.595 

-0.350 

-1.070 

-0.967 

16 

SS66444701 

82,976,515 

C/G 

-0.672 

-0.259 

0.831 

0.481 

0.639 

1.552 

0.876 

-0.143 

16 

SS66529263 

84,383,030 

G/T 

0.478 

0.539 

-0.033 

0.301 

0.066 

-0.313 

-0.037 

-0.111 

23 

SS66369851 

121,966,143 

G/T 

-0.419 

-0.108 

-0.052 

-0.025 

0.050 

-0.663 

-0.454 

-0.185 
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Fig. 2. Additive effects of selected SNPs in the real data example. 


introduced by assigning multivariate Laplace priors to regression coefficients, 
and are implemented on the basis of its hierarchical expansion which yields 
an efficient Gibbs sampler in the MCMC estimation. Although computa¬ 
tionally intensive, it outperforms the standard group lasso in the sense that 
it provides not only point estimates but also interval estimates of all pa¬ 
rameters. In addition, the Bayesian group lasso treats the regularization 
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Fig. 3. Dominant effects of selected SNPs in the real data example. 
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parameters as unknown hyperparameters and estimates them along with 
other parameters. This technique avoids choosing the tuning parameters by 
cross-validation and automatically accounts for the uncertainty in its selec¬ 
tion that affects the estimates of regression coefficients. 

In one of the most powerful but challenging areas in genetics, we incor¬ 
porated our new procedure to genome-wide association studies (GWAS) by 
testing a large number of SNPs simultaneously, particularly with n, 
based on the dynamic pattern of genetic effects on complex phenotypes 
or diseases. We first applied the new approach to /GWAS for age-specific 
changes of BMI and successfully identified several significant SNPs, some 
of which are confirmed by empirical genetic studies [Prayling (2007)]. For 
example, previous molecular studies have observed a candidate gene {FTO) 
coding alpha-ketoglutarate-dependent dioxygenase, a fat mass and obesity- 
associated protein. Our model detected SNPs ss66091573, ss66242525 and 
ss66489647 on chromosome 16 in a region of the FTO gene, suggesting the 
biological relevance of these SNPs in fat-related trait control. Our model 
also detected other SNPs in close proximity of different candidate genes; 
that is, SNP ss66397464 in peroxisome proliferator-activated receptor -7 gene 
(PPARG) on chromosome 3, SNP ss66402098 in the Wolfram syndrome 
1 gene (WFSI) on chromosome 4, SNP ss66218814 in CDK5 regulatory- 
subunit-associated protein 1-like 1 gene (CDKALl) on chromosome 6 , and 
SNP ss66288005 in potassium inwardly-rectifying channel, subfamily J, mem¬ 
ber 11 gene (KCNJll) on chromosome 12 [Frayling (2007)]. Among these 
four genes, PPARG and KGNJ were found to be associated with obesity 
[Vidal-Puig et al. (1997); Morgan et al. (2010)], while WFSI and GDKALl 
are believed to be associated with diabetes [Sandhu et al. (2007); Scott et al. 
(2007); Steinthorsdottir et al. (2007)]. Therefore, all these discoveries have 
well validated the biological relevance of the new model. 

To address challenges for the post-GWAS era, genetic association studies 
began to focus on SNPs within a set of functional candidate genes. For in¬ 
stance, Michel et al. (2010) analyzed 566 SNPs from 14 candidate genes that 
are believed to be associated with asthma. Xu and Taylor (2009) developed 
tools to recommend SNPs based on information on gene expression stud¬ 
ies, regulatory pathways and functional regions that appear to be linked 
to the disease. In their example, 1361 SNPs were recommended for a ge¬ 
netic association study on prostate cancer. These tools could be used as a 
preprocessing step for the proposed procedure in this article. Statistically, 
on the other hand, variable screening approaches [Fan and Lv (2008)] for 
longitudinal data can be developed to recommend a subset of SNPs. 

From a theoretical point of view, the proposed method can also approxi¬ 
mate varying-coefficients by nonparametric techniques other than Legendre 
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polynomials, and model the within-subject correlation by other paramet¬ 
ric or nonparametric covariance structures. Given its potential influence, 
an optimal model for longitudinal covariance structure should be chosen 
based on the nature of practical data [Zhao et al. (2005); Yap, Fan and Wu 
(2009)]. More generally, it can be easily extended to the problem where the 
number of variables in each group varies, such as the multi-factor ANOVA 
with each factor having several levels. Also, gene-gene interactions and gene- 
environment interactions can be incorporated to better decipher a detailed 
picture of the genetic architecture of a complex trait. 

Acknowledgments. Jiahan Li and Zhong Wang contributed equally to 
this work. The authors are grateful to the Editor, the Associate Editor and 
three anonymous referees for providing valuable comments that significantly 
improved the paper. 


SUPPLEMENTARY MATERIAL 

Convergence diagnostics and summary of parameter estimates (DOL 
10.1214/15-AOAS808SUPP; .pdf). We plot the potential scale reduction fac¬ 
tor (PSRF) against iterations and summarize the average estimates, stan¬ 
dard errors and mean squared errors (MSEs) of corresponding Legendre 
coefficients for the first five genetic predictors. 
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